---
title: "NVSQ_CF_DAF_code"
author: "Viviana Wu"
date: "2024-07-28"
output:
  word_document: default
  html_document: default
  pdf_document: default
---


```{r setup, include=FALSE}
library(car)
library(dplyr)
library(psych)
library(MASS)
library(sandwich)
library(stargazer)
library(mfx)
options(scipen=999, digits = 3)
```

Regression analysis
```{r}
df <- read.csv("nvsq_cf_daf_replication_data.csv")
summary(df$daf_grant5)
describe(df$daf_grant_1)

summary(m1a <- glm.nb(all_codes ~ daf_grant_1,
                     data = df))

summary(m1b <- glm.nb(all_codes ~ daf_grant_1 + 
                       TotalAssets_million + 
                       GoverningBodyVotingMembersCnt + age + SAcount,
                     data = df))

summary(m1c <- glm.nb(all_codes ~ daf_grant_1 + 
                       TotalAssets_million + 
                       GoverningBodyVotingMembersCnt + age + SAcount +
                       median_avg	+	bachelor_avg + sk_avg + 
                       pop_sum + urbanity_avg,
                     data = df))
```

Regression table
```{r}
# marginal effect
# beta = 0.2: An individual with one more year of education is predicted to have 0.2 more visits, other things equal.

m1a.m <- negbinmfx(m1a, data = df)
m1b.m <- negbinmfx(m1b, data = df)
m1c.m <- negbinmfx(m1c, data = df)

m1a$AIC <- AIC(m1a)
m1b$AIC <- AIC(m1b)
m1c$AIC <- AIC(m1c)

m1a$BIC <- BIC(m1a)
m1b$BIC <- BIC(m1b)
m1c$BIC <- BIC(m1c)

```



# Online appendices
```{r}
summary(m1 <- glm.nb(all_codes ~ daf_grant_1 + 
                       TotalAssets_million + 
                       GoverningBodyVotingMembersCnt + age + SAcount +
                       median_avg	+	bachelor_avg + sk_avg + 
                       pop_sum + urbanity_avg,
                     data = df))
vif(m1)
summary(m2 <- glm.nb(strategy_ref ~ daf_grant_1+ TotalAssets_million + 
                           GoverningBodyVotingMembersCnt + age + SAcount + 
                       median_avg	+	sk_avg + bachelor_avg +  
                       pop_sum + urbanity_avg,
                     data = df, control=glm.control(maxit=50)))

vif(m2)

summary(m3 <- glm.nb(convene_ref ~ daf_grant_1+ TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                       median_avg	+	sk_avg + bachelor_avg + 
                       pop_sum + urbanity_avg,
                     data = df))

summary(m4 <- glm.nb(knowledge_ref ~ daf_grant_1+ TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                      median_avg	+	sk_avg + bachelor_avg + pop_sum + urbanity_avg,
                     data = df))

summary(m5 <- glm.nb(capacity_ref ~ daf_grant_1 + TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                        median_avg	+	sk_avg + bachelor_avg + pop_sum + urbanity_avg,
                     data = df))

summary(m6 <- glm.nb(partner_ref ~ daf_grant_1+ TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                        median_avg	+ sk_avg + bachelor_avg + 
                       pop_sum + urbanity_avg,
                     data = df))

summary(m7 <- glm.nb(policy_ref ~ daf_grant_1 + TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                       median_avg	+ sk_avg + bachelor_avg + 
                       pop_sum + urbanity_avg,
                     data = df))

summary(m8 <- glm.nb(policy_ref ~ DAF_ind + TotalAssets_million +
                       GoverningBodyVotingMembersCnt + age + SAcount + 
                       median_avg	+ sk_avg + bachelor_avg + 
                       pop_sum + urbanity_avg,
                     data = df))


m1.m <- negbinmfx(m1, data = df)
m2.m <- negbinmfx(m2, data = df, control=glm.control(maxit=50))
m3.m <- negbinmfx(m3, data = df)
m4.m <- negbinmfx(m4, data = df)
m5.m <- negbinmfx(m5, data = df)
m6.m <- negbinmfx(m6, data = df)
m7.m <- negbinmfx(m7, data = df)

```


Summary statistics 
```{r}
summary_data2 <- round(df[c("all_codes",
                                 "strategy_ref","convene_ref","knowledge_ref",
                                 "capacity_ref","partner_ref","policy_ref", 
                                 "DAF_ind", "ScheduleD.DonorAdvisedFundsHeldCnt", 
                                 "ScheduleD.DonorAdvisedFundsVlEOYAmt",
                                 "ScheduleD.DonorAdvisedFundsGrantsAmt", 
                                 "daf_grant_1", "grants",
                                 "TotalAssets_million",
                                 "GoverningBodyVotingMembersCnt", "age", "SAcount", 
                                 "median_avg","sk_avg","bachelor_avg", 
                                 "pop_sum", "urbanity_avg")],2)

```

Summary descriptions for DAF
```{r}
# for summary stats: extract regression sample
summary(m1_daf <- glm.nb(all_codes ~ daf_grant_1 + DAF_ind + 
                           ScheduleD.DonorAdvisedFundsHeldCnt + 
                       TotalAssets_million + 
                       GoverningBodyVotingMembersCnt + age + SAcount +
                       median_avg	+	bachelor_avg + sk_avg + 
                       pop_sum + urbanity_avg,
                     data = df))

# DAF account using regression sample 
df_logic_2 <- m1_daf$model %>%
  mutate(logic2 = case_when(
    ScheduleD.DonorAdvisedFundsHeldCnt > 0 & all_codes > 0 ~ 1,
    ScheduleD.DonorAdvisedFundsHeldCnt == 0 & all_codes > 0 ~ 2,
    ScheduleD.DonorAdvisedFundsHeldCnt > 0 & all_codes == 0 ~ 3,
    ScheduleD.DonorAdvisedFundsHeldCnt == 0 & all_codes == 0 ~ 4,
    TRUE ~ 0 # default case if none of the above conditions are met
  ))
summary(factor(df_logic_2$logic2))

#  1   2   3   4 
# 1 = 392 (81.33%) 
# 2 = 19  (3.94%) 
# 3 = 64  (13.28%)
# 4 = 7 (1.45%)

```


